Electron current drive by fusion-product-excited lower hybrid drift instability 
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We present first principles simulations of the direct collisionless coupling of the free energy of 
fusion-born ions into electron current in a magnetically confined fusion plasma. These simulations 
demonstrate, for the first time, a key building block of some "alpha channelling" scenarios for toka- 
mak experiments. A fully self-consistent electromagnetic 1D3V particle-in-cell code is used to evolve 
a parallel drifting ring-beam distribution of 3MeV protons in a lOkeV thermal deuterium-electron 
plasma with realistic mass ratio. Collective instability gives rise to electromagnetic field activity 
in the lower hybrid range of frequencies. These spontaneously excited obliquely propagating waves 
undergo Landau damping on resonant electrons, drawing out an asymmetric tail in the distribution 
of electron parallel velocities, which carries a current. 



Optimal exploitation of the free energy of fusion prod- 
ucts, for example the alpha-particles born at 3.5MeV in 
reactions between thermal deuterons and tritons at 10- 
20keV, is central to achieving fusion power through mag- 
netic confinement of plasma. In the traditional frame- 
work, this energy is transferred through multiple colli- 
sions to the thermal electrons on a cumulative timescale 
~ls; electron heating by fusion alpha-particles has been 
observed in the TFTR[1] and JET [2] tokamaks. The 
electrons in turn sustain the temperature of the thermal 
ions to which they are collisionally coupled. It may be 
preferable, however, to use some of this fusion product 
free energy in alpha channelling scenarios. This term, 
coined by Fisch and Rax[3|, refers to mechanisms by 
which rapid collisionless collective instabilities (natural 
or induced) of the fusion product population could di- 
rectly benefit the plasma equilibrium, for example by 
helping sustain toroidal current . Here we report 
particle-in-cell (PIC) simulations of fusion-born protons 
in deuterium plasmas that demonstrate from first prin- 
ciples, for the first time, a key alpha channelling phe- 
nomenon for tokamak fusion plasmas. We focus on the 
collective instability of centrally born fusion products on 
banana orbits at the plasma edge, a population known 
to be responsible for observations of ion cyclotron emis- 
sion in JET[7] and TFTR[||. A fully self-consistent elec- 
tromagnetic 1D3V PIC code evolves a parallel drifting 
ring-beam distribution of 3MeV protons in a lOkeV ther- 
mal deuterium-electron plasma with realistic mass ratio. 
Collective instability gives rise to electromagnetic field 
activity in the lower hybrid range (LH) of frequencies. 
The spontaneously excited obliquely propagating waves 
undergo Landau damping on resonant electrons, draw- 
ing out an asymmetric tail in the distribution of electron 
parallel velocities, which carries a current. These sim- 
ulations thus demonstrate a key building block of some 
alpha channelling scenarios: the direct collisionless cou- 
pling of fusion product free energy into a form which can 
help sustain the basic equilibrium of the tokamak plasma. 

Spatially localized inversions of the velocity distri- 



bution of fusion-born ions can arise due to the parti- 
cle energy and pitch angle dependence of particle or- 
bits. Alpha-particles born with pitch angles just in- 
side the trapped-passing boundary generate ion cyclotron 
emission @,H|, for example. This motivates our model[9] 
of the initial fusion product velocity distribution function 
as a ring travelling anti-parallel to the magnetic field with 
distribution function f p = l/(27rv r )5(v\ \ —u)5(v± — v r ); u 
is the magnetic field aligned velocity and v r is the perpen- 
dicular velocity. The radial extent of the emitting region 
is of order a few energetic particle gyroradii in JET and 
TFTR@,[i|. 

Our simulations use physical mass ratios to ensure 
that the physically relevant instability is excited. We 
select parameter values similar to those in the edge of 
a large tokamak, subject to computational resource con- 
straints. The minority energetic fusion product protons 
at 3MeV have pitch angle 135° and their number den- 
sity n p is 1% of that of the background (Maxwellian) 
lOkeV deuterons rid- The electron number density of the 
initially quasineutral plasma is n e = 10 18 m -3 and the 
electron beta (3 e = fid = 3 x 10 -4 . The speed of the ener- 
getic protons is approximately half the Alfven speed Va 
given the applied magnetic field B = 3T. These param- 
eter values imply a total energy of the energetic protons 
~ 1.7 times that of the thermal deuterons and electrons 
combined. While this is ~ 10 times the value anticipated 
for next step fusion plasmas [lOj, it is necessary so as to 
drive the instability on an acceptable timescale compu- 
tationally. 

The collective instability of the energetic protons can 
be identified as a form of lower hybrid drift instability 
(LHDI), a topic of widespread relevance to space and 
astrophysical plasma s [Til [l2| , while LH waves are used 
for current drive [l3l-ll5| in the fusion context and are 
associated with acceleration mechanisms in ionospheric 
plasmas [l6j. The treatment here extends LHDI theory 
in several respects: the energetic ion population does not 
contribute to the plasma equilibrium, unlike most space 
and astrophysical applications, where the ion beams are 
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typically associated with currents and gradients central 
to the equilibrium; values of the key dimensionless pa- 
rameters are guided by large tokamak edge plasma con- 
ditions; and the physically correct mass ratio of electrons 
to ions is used, which was not possible in some of the clas- 
sic computational studies of LHDI. 

Our fully self-consistent, relativistic kinetic simulations 
use a particle-in-cell code epochld based on the approach 
of Ref. [I?} • Computational macroparticles represent the 
particle distribution functions in full three dimensional 
velocity space. All three components of all vector quanti- 
ties, that is, particle velocities and electromagnetic fields, 
are functions of a single spatial coordinate (x, referred to 
as the direction of variation or the simulation domain) 
and time t. Field and particle boundary conditions are 
periodic. Wavevectors are then parallel or antiparallel 
to the x direction. To focus on obliquely propagating 
modes, we set the background field at an angle = 84° 
to the x direction. The fields and bulk plasma parame- 
ters are resolved (in x) by computational cells of width 
Ax = Ad/10, where A^ is the electron Debye length. 
The simulation domain (in x) is optimized to capture 
the rather broad range of characteristic lengthscales of 
the energetic protons, background deuterons and elec- 
trons. Each of the Nq = 2048 cells corresponds to species 
gyroradii p of ~2 — 3p p , ~ 10 pd and ~ 1000p e . The simu- 
lation requires over 2 million computational particles to 
give good phase space resolution and is typically iterated 
over ~8 x 10 5 timesteps. 
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FIG. 1. Time evolution of: the total kinetic energy e% in each 
plasma species; the energies in the electric (magnetic) field 
Se (£b); the proton fluctuation energy e p . 

An overview of the energy flows in the simulation is 
shown in Figure 1. The total electric and magnetic en- 
ergy in excited fields in the simulation are obtained by 
summing the (suitably normalized) energy densities over 
all grid cells giving electric se = ^Z^Jf an d magnetic 
£b = J2(Bf ~ Bq) field energies where Bq is the applied 
magnetic field. Both the electric field energy se (light 
blue trace) and the magnetic field energy sb (magenta 
trace) rise with time. Four stages of the evolving simula- 
tion are indicated by vertical lines (i)-(iv) corresponding 
to times 8.9, 12.0, 15.2 and 17.8 lower hybrid oscillation 



periods t l h respectively. Stage (i) corresponds to the 
onset of the linear phase of field growth, which is well 
developed by stages (ii) and (iii). By stage (iv) the wave 
amplitude is approaching its saturated level. The mag- 
netic fluctuations contain nearly an order of magnitude 
less energy than the electric fluctuations, implying that 
the waves produced by the instability are largely electro- 
static. 

The total kinetic energy of each particle species is ob- 
tained by summing over all the computational particles 
in the system. From stage (iii) onwards, Figure 1 shows 
that the total kinetic energies of the protons and electrons 
approximately mirror each other, the protons losing en- 
ergy as the electrons gain energy, whilst the deuterons 
show little change. 
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Parallel Momentum [p n 1 

FIG. 2. Electron parallel momentum distribution function at 
four snapshots in time (i) — (iv) (see Figure 1) in units of ini- 
tial electron rms momentum po, e - Vertical solid black traces 
indicate the phase velocities of the fundamental modes shown 
in Figure 3, in non-relativistic electron momentum space. 

To obtain the relative gain, or loss, in the energy of 
each population we define a change in total kinetic en- 
ergy of the ensemble of particles Asi =< Si{t) > — < 
6i(t = 0) >, and a total kinetic energy of fluctuations 
€i(t) =<\sj(t)— <£i(t = 0) > | > relative to the initial 
conditions, where < ... > indicates an average over all the 
particles of species i. The time variation of the proton 
fluctuation energy e p is shown in Figure 1 . It grows from 
the start of the simulation, ultimately increasing by three 
orders of magnitude, whereas the total proton kinetic en- 
ergy declines by much less than one order of magnitude. 
This reflects the role of proton energy dispersion in the 
early phase of the instability. The electron kinetic energy 
e e rises with the electric field energy e e during the linear 
stage of the simulation, arising from electron participa- 
tion in the principally electrostatic waves excited by the 
instability. The corresponding effect for the deuterons is 
also visible, but is much less due to their higher mass. 

We now turn to the evolving distribution function of 
the electron parallel momentum p\\ shown in Figure 2 
for times (i)-(iv). From stage (iii) onwards the electron 
distribution function develops an asymmetric tail in py 
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FIG. 3. Fast Fourier transform of the E x electric field over 
the spatial domain and the time interval 10 < t/rLH < 18. 
Colour denotes field power in arbitrary units. The electron- 
deuteron cold dispersion relation (black dash trace) and the 
positions of the v x proton distribution function peaks at t = 
and time (iv) (black dash-dot traces) are overplotted. 
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reflecting net directional electron acceleration. We infer 
Landau damping of the excited waves on resonant elec- 
trons which results in the flattening of the negative tail 
of the pn distribution function. 

This is confirmed by the analysis shown in Figure 3 
where we plot the spatio-temporal fast Fourier trans- 
form of the electric field component along the simula- 
tion direction E x (uj,k), transformed from E x (x,i) for 
the time interval 10 < t/rLH < 18 ending at time (iv). 
The oblique cold-plasma normal mode of the background 
deuteron-electron plasma is marked by the black dashed 
trace. The presence of the additional energetic proton 
population modifies the nature and dispersion properties 
of the normal modes of the plasma and, through reso- 
nance, couple energy to these modified modes. The Fig- 
ure shows peaks in intensity at (lj, k) ~ (Qlulh, ±2u pe /c) 
which indicate resonance of the proton population with 
a modified deuteron-electron normal mode branch that 
lies on the surface between the lower extraordinary wave 
and the whistler wave in LJ,k±,k\\ space. The locations 
in lj, k space where coupling takes place correspond to 
where the velocities of the peaks of the proton distribu- 
tion function in v x approximately match the phase ve- 
locities of the normal mode. As the proton distribution 
function evolves, the positions in v x of the peaks in the 
distribution function move. These velocities at t = and 
at time (iv) are plotted on Figure 3 as black dash-dot 
traces. We note that there are harmonics of the dominant 
backward propagating (in negative x direction) mode at 
(lj, k) ~ (6u LH , -2u pe /c). 

On Figure 2 we indicate with vertical lines the two 
electron parallel momenta resonant with the dominant 
backward and forward electromagnetic structures shown 
in Figure 3. Taken together, Figures 2 and 3 indicate that 
the electrons are principally accelerated by the dominant 
wave at negative phase velocity, which is in turn excited 
by the backward-drifting fusion products. 
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FIG. 4. Spatio-temporal plots of the backward E~ (x, t) (left) 
and forward E x (x,t) (right) travelling wave components of 
E x (x,t). Colour indicates normalized field amplitude. 

To capture the coherent oscillatory features of the 
dominant excited fields, and hence their phase resonant 
characteristics, we decompose the electric field E x (x,t) 
into component waves travelling with positive and nega- 
tive phase velocities. These are E+(x,t) = IFT[E x (lj> 
0,fc>0) + E x (lj <0,k<0)] and E~(x,t) = IFT[E x (lj< 
0,k> 0)+E x (lj > 0, k < 0)] , where IFT denotes the inverse 
Fourier transform, which we plot in Figure 4. The sum 
of E+(x,t) and E~(x,t) recreates the original field. In 
Figure 4 , the amplitude of the backward travelling wave 
is approximately an order of magnitude greater than the 
forward travelling wave and is dominated by a single long 
wavelength component. The parallel phase velocity of 
this backward travelling wave is shown by the vertical 
line on the left of Figure 2. 

The predominantly drift character of the instability is 
directly seen in the fully resolved velocity space of the 
PIC simulation. The velocity of each particle can be 
expressed in terms of a field aligned component vn, a 
component aligned with the simulation domain (and k) 
v X: and gyrophase arctan(v±^/v±^)- We test for veloc- 
ity space patterns in these coordinates at snapshot (hi) 
when the wavefields are well established. We select pro- 
tons from a narrow region in configuration space Sx that 
is smaller in extent that the wavelength of the dominant 
wavemode (Sx/X = 0.1). The coordinates v\\ : v x , and 
the gyrophase for each particle are plotted in panel (a) 
of Figure 5, where gyrophase is along the abscissa, vu is 
along the ordinate and colour indicates v x . We see asym- 
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FIG. 5. Panel (a): Snapshot at time (iii) of the energetic 
proton velocity space. Velocity space coordinates are v\\ (ab- 
scissa) , gyrophase (ordinate) and v x (colour) . Panel (b) : Nor- 
malized wave amplitude seen by protons at resonance with the 
dominant wave plotted as a function of phase space (see text). 
Panel (c): As in panel (b) for a sum of dominant wavemode 
and the counter-propagating damped mode. 

metry in the oscillatory pattern in v\\ which is a function 
of v x : a characteristic of drift instability rather than gy- 
roresonance. As such, these oscillatory perturbations in 
vn should vary as the wave amplitude at the point of 
resonance v x = ou/k, and in panel (b) of Figure 5 we con- 
firm that there is indeed resonance with the dominant 
LH wave identified above. Panel 5(b) plots on the ordi- 
nate the normalized wave amplitude § , that is, the sine 
of the wave phase (kxR — gj£r) at the location xr and 
time tR of unperturbed test particles that satisfy the con- 
dition for resonance with the dominant wave (gj, fc). The 
point of resonance (xr^r) is obtained for test particles 
which initially are distributed uniformly in gyrophase <fi 
and which follow unperturbed cycloidal orbits. To recon- 
struct a single snapshot in x and t (to compare with the 



simulation snapshot of panel (a)) these test particle orbits 
are then advanced in x, v x and to a single point (#1, ti). 
We plot in panel (b), for a snapshot in space and time 
(xi,ti), the normalized wave amplitude experienced by 
the particle at the point of resonance <o{xr^r) against 
particle gyrophase <j>(t\) with v x ((j)(ti)), represented in 
colour. We see that panel 5(b) closely tracks the large 
scale oscillations seen in the PIC code velocity space of 
panel 5(a). We refine this idea in panel 5(c) where we 
plot the effective wave amplitude arising from both the 
forward and backward waves identified in Figure 4; the 
weaker second wave can be seen to introduce fine scale 
oscillations visible in panel 5(a). 

These first principles simulations demonstrate, for the 
first time, key physical elements of alpha channelling sce- 
narios for future tokamak plasmas: the excitation of LH 
waves by a fusion product population, whose functional 
form and parameters are aligned with prior observations 
on JET and TFTR, combined with the subsequent Lan- 
dau damping of the excited waves on resonant electrons, 
drawing out an asymmetric tail in their parallel veloc- 
ity distribution, which carries a current. These simula- 
tions also contribute to the question of what instabili- 
ties may arise - transiently, locally, or otherwise - during 
the initiation and propagation of fusion burn in tokamak 
plasmas. Furthermore they deepen our understanding of 
LH drift instabilities, which may be widespread in space 
and astrophysical plasmas, by extending their study into 
parameter ranges that approximate (subject to compu- 
tational resource constraints) edge plasma conditions in 
large tokamaks. 
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